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Abstract. The quantum dynamics of interacting bosons in a one-dimensional system 
is investigated numerically. We consider dissipative and conservative two-particle 
interactions, and integrate the master equation describing the system dynamics via a 
timc-cvolving block-decimation (TEBD) algorithm. Our numerical simulations directly 
apply to stationary-light polaritons in systems where atoms and photons are confined 
to the hollow core of a photonic crystal fibre. We show that a two-particle loss term 
can drive an initially uncorrelated state into a regime where correlations effectively 
inhibit the dissipation of particles. The correlations induced by two-particle losses 
are compared with those generated by an elastic repulsion. For the considered time 
range, we find a similar behaviour in local density-density correlations but differences 
in non-local correlations. 
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1. Introduction 

The preparation of quantum systems in entangled many-body states is an indispensable 
resource for quantum information science [1] and allows one to access the rich physics of 
strongly correlated many-body systems [2]. Despite its reputation as the main adversary 
of quantum information schemes, dissipation was recognised recently as a versatile tool 
for the generation of specific many-body states [3], HI [5] . 

The preparation of correlated quantum states via dissipation requires that the 
target state itself is metastable. An intriguing experiment [6] with cold molecules 
showed that a dissipative two-particle interaction between bosons in ID gives rise to 
correlations that inhibit the loss of particles. The formal analysis of this effect is based 
on a generalised Lieb-Liniger model [3 [3 [9] and suggests that the two-particle loss 
term effectively results in a repulsion such that inelastic collisions between particles are 
avoided. In the limit of strong interactions [8], the two-particle losses give rise to a 
Tonks-Girardeau gas [10] where bosons behave with respect to many observables as if 
they were fermions. 

Recently, it was shown that ID systems of interacting bosons can be realised 
via dark-state polaritons that arise in light-matter interactions under conditions of 
electromagnetically induced transparency [TT| [T2| [T3| [HI [15]. A possible experimental 
realisation of the polariton schemes is the setup described in [16], |T7], where photons 
and atoms are simultaneously confined to the hollow core of a photonic-crystal fibre, 
see Fig. [H Another approach comprises the experimental setup in [18], where 
atoms are trapped near the surface of an optical nanofibre using evanescent waves. 
The nature of the polariton-polariton interaction in these approaches can be tuned 
by external parameters and can be elastic [I3] or dissipative [HI |T5], and the 
interaction strength is maximal for a purely dissipative interaction [T3]. Polariton 
schemes [13], IHl [TS] [191 ED] 121] are extremely promising candidates for the realisation 
of strongly correlated photon states. 
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Here we consider bosons in ID and investigate the quantum dynamics after 
a sudden switch-on of elastic and inelastic interactions via time-dependent density 
matrix renormalisation group techniques. The master equation of the system and the 
uncorrelated initial state are chosen such that our simulations directly apply to polariton 
systems [I3l [HI [15]. All details of our model are discussed in Sec. [21 and numerical 
results are presented in Sec. [3] In particular, we are interested in the preparation of a 
correlated, non-decaying state via dissipative two-particle interactions that we discuss 
in Sec. 13.11 Non-local correlations are investigated in Sec. 13.21 where we show that 
dissipation and repulsion give rise to different results. Finally, a summary is provided 
in Sec. [H 



2. Model 



We consider bosons of mass m in a one-dimensional system of length L that experience a 
two-particle contact interaction that can be either conservative, dissipative or a mixture 
of both. First we describe the full master equation governing the time evolution of 
the system in Sec. 12.11 A physical realisation of this model via interacting dark-state 
polaritons [HI [15] is introduced in Sec. 12.21 Finally, we discuss the discretised model 
of the master equation in Sec. 12.11 that allows us to perform numerical simulations in 
Sec. [231 



2.1. Master equation 

The dynamics of bosons in our ID system is described by the master equation 

hdtQ = -iHcsQ + igHlg +Ig + hCriQ, (1) 
where Hes is a non-hermitian Hamiltonian, 

Here ip{z) is the field operator obeying bosonic commutation relations 

[i;iz),i;^{z')] = 5{z-z'), (3) 

and the first term in Eq. (|2]) describes the kinetic energy. The parameter g is the 
complex coupling constant and its real part gives rise to a hermitian contribution to 
Hes in Eq. ([2]) that accounts for elastic two-particle collisions. On the other hand, the 
imaginary part of g together with the term 

Ig = -lm{g) [ dzip^gip^^ (4) 



in Eq. constitutes a two-particle loss term that can be written in Lindblad form [22] 
as lm{g/2)V[^lj% where V[X] = J^dz{X^Xg + gX^X - 2XgX^). The last term in 
Eq. ([T]) is defined as 

Cug = / dzV[d,^] (5) 
^ Jo 
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and causes diffusion of the field if). Throughout this paper we choose 

D = \h/m\/10 (6) 

such that the losses associated with the diffusion term are small on a timescale set by 
the kinetic energy term in Eq. ([2]). 

The master equation ([1]) without the last term Cdq defined in Eq. ([5]) can 
be identified with the dissipative Lieb-Liniger model [8]. For a fixed number of 
particles A^ph, essential features of the system are described by the complex Lieb-Liniger 
parameter 

_ mg 

^ - Wn-Jl- 

The spectrum of the effective Hamiltonian in Eq. ([2]) was analysed in [S] via the Bethe- 
Ansatz. It was found that the ground state of the so-called gaseous states (all momenta 
converge to finite values for |G| — >■ c>o) approaches the Tonks- Girardeau gas in the limit 
|G| — )■ oo. This state is metastable since in the Tonks- Girardeau gas limit two particles 
never occupy the same position in space, and hence the two-particle contact interaction 
does not contribute to the dissipation of particles. 

This feature can be formally derived from the master equation ([1]). The two-particle 
loss term lm{g/2)V[ip'^] in Eq. ([T]) causes a time dependence of the density of particles 
according to 

dt{h{z)) = ^lm{~g)g^'\z,z){n{z)f, (8) 
where n(2) = ip'^{z)ilj{z) is the particle density operator and 

g yz,z)— \\/-/ /\\ \^) 

is the second order correlation function. Note that the kinetic energy and the diffusion 
term Cdq in Eq. ([1]) will, in general, induce an additional time dependence of {n{z)). The 
kinetic energy term conserves the total number of particles but gives rise to polariton 
fluxes from and to position z that would redistribute the polariton density. In regions 
where the system and the particle density are homogeneous, these fluxes can safely 
be neglected. The impact of the diffusion term CbQ is negligible for the considered 
time scales since we chose a small value for the parameter D in Eq. ([6]). In the 
strongly correlated regime and for a homogeneous system, the gaseous ground state 
of iJcff obeys [S] 

^/(2)(^,^) = (l-l/Ar2j47rV(3|Gp), (10) 

and hence g^'^\z, z) -C 1 for |G| sufficiently large. The combination of Eqs. (IHl) and (fTOl) 
implies that a purely dissipative interaction term supports metastable states. 
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field n+ 

Figure 1. (a) Considered setup of A^a atoms confined to an interaction volume of 
length L and transverse area A. Q,± are the Rabi frequencies of the classical control 
fields, and E± are the quantum probe fields, (b) Atomic level scheme. 7ij is the full 
decay rate on the \i) o \j) transition, 5 and A label the detuning of the probe fields 
with states |3) and |4), respectively, and e is the two-photon detuning. 

2.2. Physical realisation 

An example for a physical system that is described by the master equation ([1]) is 
shown in Fig. [I^a), where photons and atoms are simultaneously confined to the hollow 
core of a photonic-crystal fibre. Since the light-guiding core of the optical fibre is of 
the same order of magnitude as the optical wavelength, the fibre represents a one- 
dimensional waveguide for the optical fields. The level scheme of each atom inside the 
fibre is shown in Fig. [U^b). The dynamics of the system can be described in terms of 
dark-state polaritons [23] representing bosonic quasi-particles that arise in light-matter 
interactions under conditions of electromagnetically induced transparency (KIT) [2^ . 
The two counter-propagating control fields VL± [see Fig. [T^a)] of equal intensity give 
rise to stationary light |25l |26] where the polaritons experience a quadratic dispersion 
relation like bosons in free space. Moreover, the coupling of the probe fields to the 
transition |2) ^ |4) induces a two-particle contact interaction between the polaritons 
that can be conservative |TT1 [13], dissipative or a mixture of both [T^ [T5] . It can be 
shown [m [15] that the dark-state polariton dynamics can be described by the master 
equation ([1]), where the effective mass m and the complex coupling constant g depend 
on the detunings 5, A [see Fig. [H^b)] and the intensity of the control fields. In particular, 
we point out that a proper adjustment of the detuning 5 allows one to fulfil Eq. (|6]), and 
the choice of the detuning A determines the relative contribution of elastic and inelastic 
processes to the two-particle interaction. 

In an experiment, the evolution of the polaritons under Eq. ([1]) has to be preceded 
by a loading process that could be realised as follows [13]. Initially the probe field 
modes are empty and only the control field is switched on. Then a probe pulse 
copropagating with enters the fibre under conditions of electromagnetically induced 
transparency. Here we assume that the transition |4) -H- |2) is far off-resonant during 
the loading process such that its influence is negligible. The dynamics of the probe fleld 
inside the medium can be described concisely in terms of dark-state polaritons [23]. As 
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compared to the pulse propagating in vacuum, the polariton pulse inside the medium is 
spatially compressed and travels with a reduced group velocity Vg oc fi^. Once the pulse 
is entirely inside the fibre, the group velocity of the polariton pulse is reduced to zero by 
an adiabatic switch-off of the control field The latter process maps the polaritons 
into a stationary spin coherence, and all probe field modes are empty. Eventually both 
control fields are switched on adiabatically and with equal intensity, VL^ = Qq. The time 
evolution of the corresponding dark-state polaritons is then determined by the master 
equation ([1]) [T^ [T5] . Note that the adiabatic switch-on can be much faster than all 
typical timescales in the master equation ([1]). At this stage, the detuning A of the probe 
field with the |4) ^ |2) transition can be adjusted by a proper choice of the frequency 
of the counterpropagating control fields. This frequency can be different from the one 
used for the single control field during the loading process. 

The correlations that built up under the evolution of the master equation ([T]) can be 
measured if the intensity of the control field f2_ is reduced. The latter process converts 
the stationary polaritons into a propagating pulse with controllable group velocity. This 
procedure maps spatial correlations of the trapped pulse into temporal correlations of 
the output light that can be detected via standard quantum optical techniques. 



2.3. Discretised model 

The dynamics of the system is studied via the time-evolving block decimation 
algorithm [271 EHl EH] for density operators. The discretised version [30] of Eq. ([T]) 
is given by a generalised Bose-Hubbard model augmented by loss terms, 

g = -^[Hi + H2,g]+Cig + C20. (11) 

Here Hi and H2 correspond to the kinetic energy term and the conservative contribution 
to the two-particle interaction, respectively, 

Hi=2JY^ ajdi - J ^(aja!+i + fli+iSj), ^2 = ^ ^J^^^' '^^^^ 
; I I 

and J and U are defined as 

The parameter Az is the lattice constant of the discretisation grid. The terms CiQ 
and £2^ in Eq- (tl]) are the discretised versions of the diffusion term in Eq. ([5]) and the 
dissipative contribution to the two-particle interaction, respectively, 

CiQ = - Tl '^{alaiQ + gajai - 2aiQa\) (14) 
I 
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Figure 2. Spatial density distribution {n{z)) of the initial state in Eq. (|17p . 



^--29 = -— 2_^{a\ aiQ + gal ai 



and 



T,=D- 



1 



-lm{g/h) — . 



(15) 



(16) 



The discretised version Eq. ( fTTl) is a good approximation of the continuous master 
equation ([1]) if the smallest wavelength Amin involved is much larger than the grid spacing 
Az. Initially Amin is determined by the momentum components of the initial state, but 
the evolution under Eq. ( fTTl) will change the momentum distribution of the sample. 
In particular, the suddenly switched-on two-particle repulsion term leads (among other 
processes) to a population of higher momentum states. It was found [3l] that the 
lattice model is a good approximation provided that U/J<^1. If the latter inequality 
is violated, lattice artefacts occur. We find that these effects do not occur for a purely 
dissipative two-particle interaction {U = 0) even if hV2/ J ~ 1. 

For all numerical simulations we consider a grid with A'^^ = 500 sites, and the 
physical state space at each site I is spanned by the four Fock states {|0);, |2);, |3)i} 
allowing for a maximum occupation of three particles at each site. We verified the 
convergence of our algorithm by varying the bond dimension x of matrix product 
state. In our simulations, the truncation errors Cti-unc are bounded by Ctrunc < 10~^ for 
X = 100. 

The choice of the initial state for the numerical integration is motivated by the 
loading process of an optical pulse into the fibre as described in Sec. 12.21 We consider an 
initial state with a mean number of (A^)o = 5 particles and a spatial density distribution 
as shown in Fig. [2l Here we neglect distortions and correlations due to the optical 
nonlinearity that may build up during the loading process. The initial shape of the 
polariton pulse is then determined by the slowly varying envelope of the probe pulse 
entering the system [23]. Furthermore, we assume that the initial polariton pulse is 
derived from laser light and thus described by a product of coherent states. In the 
discretised model we thus model the initial state by a coherent state at each site, 

|V^)o = nexp(-|Q|V2)exp(Qa;)|0), (17) 
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Figure 3. Temporal evolution of ^^^^(zo, zq) at the centre zq = L/2 of the sample for 
(a) \G\ = 1, (b) \G\ ^ 10, (c) \G\ = 20 and (d) \G\ = 100. The timescale Tc is defined in 
Eq. ()19|) , and the solid lines correspond to a purely dissipative two-particle interaction, 
Re(p) = 0. In (a) and (b), the dashed line shows the result for an interaction that is 
predominantly repulsive [|Re(^)/Im(^)| = 10]. 



where J2i = 5 and |qP = Az{n{zi)) with zi = lAz. We point out that |'0)o in 
Eq. (fT7|) is an uncorrelated state with g^'^^z, z') = 1 for all values of z, z' in the sample. 

3. Numerical results 

Next we describe the results of the numerical integration of the master equation ( ITT]) . 
Starting from the initial state in Eq. (ITTl) . the interaction is suddenly switched on at 
t = 0. We consider various scenarios where the interaction is purely dissipative or 
predominantly conservative. In Sec. 13.11 we discuss whether two-particle losses alone 
are able to drive the system into a correlated regime where dissipation of particles 
is suppressed. As discussed in Sec. 12. the inhibition of two-particle losses is linked 
to local correlations in the sample. The second part Sec. 13.21 is concerned with non- 
local correlations. In particular, we compare non-local correlations created by a purely 
dissipative and a predominantly conservative interaction. 

3.1. Inhibition of two-particle losses 

The two-particle decay term in Eq. ([1]) in general leads to a loss of particles from the 
system. The rate at which the particle density diminishes is determined by Eq. (jH]) and 
depends crucially on the local correlations g^'^\z,z) of the system. Figure [3] shows the 
temporal evolution of g^'^\zQ, zq) at the centre zq = L/2 of the sample for different values 
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of the Lieb-Liniger parameter |G| according to the numerical integration of Eq. fllip . 
The values of IGj correspond to the mean particle number (iV)o = 5 of the initial state, 
see Eq. ( ITTl) . All solid lined lines in Fig. [3] belong to a purely dissipative two-particle 
interaction, Re(^) = 0. On the contrary, the dashed lines in Fig. [3] (a) and (b) correspond 
to a predominantly conservative interaction |Re(^)/Im(^)| = 10. 

In the case of a purely conservative interaction (Im(^) = 0), it was shown [31] that 
local quantities reach a stationary state on a timescale Tioc = ^/ {gp)i where p = N^^/L is 
the density of particles. Note that the Lieb-Liniger gas does not thermalize globally |32j . 
Here we generalise this definition to the case of complex g and define 

h m 1 
\g\p h \G\p^ 

Due to the enhanced numerical complexity of TEBD for density operators and numerical 
limitations, our simulations cannot cover the full time range Tioc. In the following we 
scale time in the presentation of numerical results with 

^ ^^loc- (19) 



The comparison between the dashed and solid lines in Figs. |3] (a) and (b) shows that 
dissipative and repulsive two-particle interactions are equally effective in the repulsion 
of particles on a timescale proportional to the inverse interaction strength |^|. Note that 
for the polariton system discussed in Sec. 12.21 \g\ is at least one magnitude larger for a 
purely dissipative interaction [TU [15] as compared to the conservative case. It follows 
that small values of g^'^\zQ, zq) are reached much faster in absolute time for a dissipative 
two-particle interaction. 

We emphasise that the minimal values for g^'^\zQ, zq) in Fig. [3] will decrease further 
in time until their equilibrium values are reached. An estimate of these equilibrium 
values can be obtained via the work presented in [31]. Here the authors performed 
TEBD simulations for the unitary time evolution of bosons in ID with two-particle 
repulsion. In particular, this system could be described by a pure quantum state which 
reduces the numerical complexity as compared to our simulations of the full density 
operator. Therefore, the simulations in [31] could be carried out up to times t = Tioc- 
Given that the initial dynamics of g^^^'^ZQ, zq) for two-particle losses and for conservative 
interactions is very similar, one may suppose that the equilibrium values of g^'^\zo, zq) 
will also be almost the same. For this case we estimate that g^'^\zQ, zq) for our system 
should lie within the range [0.1, 0.2] for |G| = 10 and |G| = 20 at t = Tioc- 

Next we determine the temporal evolution of the particle density at the centre of 
the cloud via Eq. ^ and the numerical results for the temporal evolution of g^'^\zo, zq). 
The result is represented by the solid lines in Fig. H] and in good agreement with the 
direct evaluation of the particle density {n{z)) via the numerical integration of the master 
equation (fTTl) (dotted lines). The small deviations are due to the diffusion term in Eq. ([5]) 
that results in particle losses and that were omitted in Eq. ([8]). Note that all curves in 
Fig. m correspond to a purely dissipative two-particle interaction, Re(^) = 0. In order to 
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Figure 4. Temporal evolution of {n{z)) at the centre zo = L/2 of the sample for (a) 
|G| = 1, (b) |G| = 10, (c) |G| = 20, (d) |G| = 100. In (a)-(d), the solid line corresponds 
to a numerical integration of Eq. (jH]) with g^'^'>(zo, zq) taken from the TEBD simulation 
for a purely dissipative two-particle interaction, Re(g) = 0. The dotted lines indicate 
the results of the TEBD simulation for {n(z)). The dashed line shows the evolution of 
{n{z)) according to Eq. ([8]) with g^^'>{zo, zq) = 1 for all times. 



visualise the slowdown of two-particle losses due to the decrease of g^'^\zo, zq) in time, 
the dashed line in Fig.|4]shows the temporal evolution of {h{z)) according to Eq. ([8]) with 
g^'^\zQ,ZQ) = 1 for all times. The slowdown of two-particle losses is most pronounced 
for |G| = 100 where a quasi-stationary regime is reached with g^'^^zo, Zq) close to zero. 
We find that the decrease of the mean number of particles {N)/{N)o follows the same 
curves as the particle density at the centre of the cloud. It follows that more than 
90% of the initially present particles are left at t/Tc = 4 for |G| = 100. Figures |4] (b) 
and (c) correspond to \G\ = 10 and \G\ = 20, respectively, and demonstrate a clear 
slowdown of two-particle losses. On the basis of the estimates for the equilibrium values 
of g^'^\zo, Zq) obtained via the results presented in |3l], we suppose that 50% (60%) of 
the initial particles are still present at t = Tioc for |G| = 10 (|G| = 20). On the contrary, 
g'''^\zo, Zq) will not drop to values close to zero for |G| = 1 and hence we expect that a 
significant number of particles is lost at t = Tioc. 

3.2. Non-local correlations 

The spatial dependence of g^'^\zQ, z) as a function of z (zq = L/2 is the centre of the 
sample) is depicted in Fig. [51 Several snapshots in time are shown indicating that non- 
local correlations exhibit oscillations that propagate outwards. The oscillations for the 
longest evolution times in Fig. |5] exhibit several minima that lie below unity. It follows 
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Figure 5. Spatial dependence of (^'^^'(zqjz) for a purely dissipative two-particle 
interaction [Re(g) = 0]. The parameters correspond to (a) \G\ ~ 20 and (b) \G\ ~ 100, 
and ZQ — Ljl is the centre of the sample. 




Figure 6. Spatial dependence of g'^'(zo,-z) for (a) \G\ — 1, t/rc — 0.056 and (b) 
\G\ ~ 10, t/rc = 0.77. zq ~ L/2 is the centre of the sample. The solid lines correspond 
to a purely dissipative two-particle interaction [Re(g) = 0], and the dashed lines show 
the result for an interaction that is predominantly repulsive [|Re(g)/Im(g)| = 10]. 



that a dissipative two-particle interaction alone to some extend gives rise to spatial order 
of the particles. Note that non-local correlations reach a quasi-stationary regime on a 
timescale that is much larger than Tioc- In particular, we point out that the correlations 
shown in Fig. [6] are not Friedel oscillations |2] that characterise the Fermionic nature of 
particles in the Tonks-Girardeau regime. 

In Sec. 13.11 we concluded that repulsion and dissipation are equally effective in the 
creation of small values of g^'^\z,z). On the other hand, dissipative and conservative 
two-particle interactions give rise to different non-local correlations as shown in Fig. [61 
The solid and dashed lines correspond to a purely dissipative and a mostly conservative 
two-particle interaction, respectively. It can be seen that the non-local correlations 
propagate at the same velocity in space in the dissipative and conservative case, but the 
oscillations are more pronounced for a repulsive two-particle interaction. In addition, 
two-particle losses give rise to a broader dip in g^'^\zo, z) as a function of z as compared 
to the conservative case. Note that the oscillation amplitudes in g^'^\zQ^z) are also 
damped by the diffusion term in Eq. ([5]). 
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4. Discussion and summary 

In this paper we investigated the quantum dynamics of bosons in ID after a sudden 
switch-on of two-particle losses. The master equation and the uncorrelated initial 
state are chosen such that our simulations directly apply to stationary-light polariton 
systems [13], [lH [15] . Apart from a (small) diffusion term, the considered master equation 
coincides with the dissipative Lieb-Liniger model [8]. In particular, we investigated 
the dissipation-induced preparation of correlations that effectively inhibit two-particle 
losses. We found that a metastable regime where losses are substantially suppressed can 
be prepared for large values of the Lieb-Liniger parameter \G\ ^ 1. Formally, it was 
shown [71 [8] that the two-particle loss term effectively results in a repulsion between 
particles such that they never occupy the same position in space. Therefore, dissipation 
via the two-particle contact interaction is inhibited. A physical explanation for this 
counterintuitive result can be given for a discrete lattice model [6l [7|. There it can be 
regarded as a manifestation of the quantum Zeno effect, where the two-particle losses 
play the role of a continuous measurement that always projects the system onto states 
with less than two particles in each lattice site before higher particle numbers can build 
up. Furthermore, we compared the correlations induced by dissipative and repulsive 
interactions, respectively. We found that dissipation and repulsion are equally efficient 
for the generation of local correlations g^'^\z,z) < 1 on a timescale proportional to 
the inverse interaction strength. On the other hand, non-local correlations induced by 
dissipation show different features than those created by two-particle repulsion. As 
compared to an elastic repulsion between the particles, we find that two-particle losses 
give rise to a broader dip in g^'^^ZQ, z) as a function of z [zq is the centre of the sample). 
Both repulsion and dissipation give rise to spatial oscillations in g^'^\zQ^ z) that spread in 
time and indicate to some extend a spatial ordering. We find that these oscillations are 
more pronounced for a repulsion between the particles as compared to the dissipative 
case. 
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